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Abstract. A theoretical analysis of the scattering of plane- wave atomic excitations 
in disordered solids has been made in terms of the spectral densities. Hybridization 
between transverse and longitudinal waves of approximately the same frequency is 
demonstrated. The analytic results agree well with the results obtained from computer 
simulation for a toy linear zig-zag chain model and a model of vitreous silica constructed 
by molecular dynamics. 



I INTRODUCTION 

Propagation of classical plane waves in random scattering media has attracted a 
lot of theoretical and experimental attention in recent years [1,2]. The phenomena 
of wide interest include weak localization [1,2], Anderson localization [3], phonon 
localization [4-6] and related behaviour of plane waves in the Ioffe-Regel crossover 
region [7] separating weakly and strongly scattered waves [8,9]. 

Vibrational plane waves propagate well in disordered structures in the long- 
wavelength limit, k<l, with k the wavevector and a a measure of the microscopic 
scale of the structure (being of the order of interatomic distances) , when the atomic 
structure behaves as an elastic continuum and disorder on the microscopic level is 
not of great importance. The situation is changed with decreasing wavelength and, 
on the microscopic scale, disorder becomes important and the wavevector is no 
longer a good quantum value [10-17]. 

In the investigation of the propagation of the plane-wave vibrational excita- 
tions in disordered structures, different decay channels have been suggested to 
explain their attenuation: (i) disorder-induced channels [18-22], (ii) anharmonic 
channels [20] and (iii) channels involving two-level systems [23-26]. The anhar- 
monic channels are strongly enhanced with increasing temperature, particularly 
at temperatures comparable with the glass-transition temperature, T g ~ 10 3 K. In 



contrast, scattering by two- level systems can be important at low temperatures 
[25], T < Ttls ~ 10 - 100K [24,27]. In the intermediate temperature range, 
Ttls < T <C T g , which is considered below, the scattering processes involving two- 
level systems are suppressed and the atomic dynamics are usually harmonic [21,28], 
meaning that disorder-induced channels play the most important role in the decay 
mechanism of plane-wave excitations. 

If the harmonic approximation is valid for describing the atomic dynamics, then 
a normal-mode analysis can be used for the problem under consideration. The 
normal modes can be found either analytically or numerically. A general theory 
of atomic vibrations in disordered structures, which in principle should result in 
normal modes, has been mainly developed for particular simple model structures 
[18,29-34], which can hardly describe quantitatively the situation in real struc- 
tures. Therefore, a numerical approach could be very useful in the calculation of 
normal modes, e.g. by direct diagonalization of the dynamical matrix which can 
be available, e.g. from molecular dynamics simulations. 

The main questions we address in this paper are: how are plane-wave vibrational 
excitations scattered by disorder and what are the characteristics of the final state 
after scattering? 

Our approach to the problem is based on combining analytical and numerical 
techniques. First, we create a few structural models of a disordered atomic material: 
(i) realistic models of v-SiC>2 using molecular dynamics and (ii) toy models of a 
linear zig-zag chain. Then all eigenmodes and eigenfrequencies (in the harmonic 
approximation) are found numerically. These characteristics fully determine the 
dynamical response of the system (and final state at t — > oo) to any external 
excitations, including the plane- wave excitations of present interest. The final 
state after scattering was investigated then in momentum space (analytically and 
numerically) . 

II FORMALISM 

The time evolution of any vibrational excitation is fully determined in the har- 
monic approximation by the eigenmodes and eigenfrequencies of the system. In- 
deed, the initial excitation can be expanded in the eigenmodes, the time dependence 
of which is known. The coefficients in such an expansion are defined by the shape of 
the initial vibrational excitation. Here we consider only plane-wave external initial 
excitations, mainly because exactly such excitations are generated in a system by 
inelastic neutron, light and electron scattering [35]. 

In amorphous materials, because of disorder, the eigenmodes are not plane waves 
even in the long- wavelength limit. Therefore, an initial plane wave, when expanded 
over eigenmodes, contains different eigenmodes characterized by different weights 
in this expansion. The eigenmodes participating in the expansion evolve differently 
with time, so that the propagating excitation becomes different in shape compared 
with the initial one. If we expand the vibrational state in plane waves after a certain 



time, then this expansion contains not only the initial plane-wave component but 
also other plane waves characterized by different wavevectors. This means that the 
initial plane wave is scattered by the structure into a different final state. Our aim 
here is to study the properties of the final state after decay, for different wavevectors 
of an initial plane wave. 

Let us consider an external wave excitation introduced in the system, u(t), which 
at the initial moment of time is an ideal plane wave, w^, characterized by the 
wavevector k, unit polarization vector n and initial phase O : 

u(t = 0) = w k = Ancos[k-r + o ] , (1) 

where u is a 3iV-dimensional displacement vector, the i-th component of which 
describes the displacement of atom i from its equilibrium position at r ; , A is the 
normalization constant defined below and the wavevector index k includes also the 
polarization index n. In our analytical treatment, we assume that eigenmodes and 
eigenfrequencies are known, e.g. from numerical simulations. The initial displace- 
ment vector, Eq. (1), each atomic component of which is multiplied by the mass 
factor rrii = MjiV/ J2i Mi (Mi stands for the mass of atom i) can be expanded in 
eigenmodes as: 

3 JV 

u(0) = S^ k e , /Vm, (2) 

where the symbolic script e 1 / ' \pm means that each i-th component of vector 
is divided by the factor y/ml. The coefficients a k in expansion (2), the squares of 
which are the spectral densities of the system, are defined by the following equation: 

TV 

a J k = (e ] ■ v^Uk) = V™i e i ' w k,i • (3) 
i=i 

The spectral-density coefficients, Eq. (3), fully determine the dynamical response 
of the system to plane-wave excitation. Indeed, at any moment of time t, the dis- 
placement vector of the propagating excitation can be represented via eigenmodes 
developing in time as: 

m _■ e j 

u(0 = J2^i^ cosuj J t ■ ( 4 ) 



For the sake of simplicity and without loss of generality (as shown below), we 
consider the initial excitation to be a standing wave, i.e. u(0) = 0, leading to the 
absence of terms proportional to sincc^-t in expression (4). It is convenient for the 
initial vector y / mu k (0) to be normalized to unity, so that 

3iV 

EKI 2 = i, (5) 



and the normalization constant in Eq. (1) is A 2 = [J2i ^j|u k (0)| 2 ] . 

An ideal initial plane wave (1) scatters with time to different plane waves. In 
order to calculate the weights of different plane- wave components in the propagating 
excitation, we expand the displacement vector u(t) in plane waves: 

u(t) = EM*)> (e) 

k' 

where the sum is taken over all wavevectors k' (allowed by the periodic simulation 
box in the case of a finite model) and all polarizations (two transverse and one 
longitudinal for each wavevector), and the waves u k /(£) are defined as: 

u k '(0 = a*? (t) An' cos (k' • r + k >(t)) . (7) 

The same normalization as in Eq. (1) is used here. In order to find the time 
dependence of the amplitude a k /(t) and phase (fik>(t), it is convenient to rewrite 
Eq. (7) in the following form: 

u k '(t) = ak',c(*)wk',c + a k ', s (i)w k / iS , (8) 

where 



w k' c — Ah' cos k' • r and w k ' s = Ah' sin k' • r 



(9) 



so that 



« k 'W = (4',cW + 4,sW) 1/2 , (io) 

k '(*) = Arctan[a k ', B (*)/atf >c (f)] . (11) 

The coefficients a k / ;C ( s )(t) before the cos- (sin-) like components in Eq. (8) can be 
found by multiplying both sides of this equation by Wk' jC ( s ) an d using Eq. (3), so 
that 

a k ',c(s)(0 = Yl «k«k',c(s) cosujt/ (w k , iC(s) ) . (12) 
j 

where 

N 

( W k',c(s)) =EK W k',c(s))i| 2 • (13) 
i=l 

The spectral-density coefficients a k , c ^ s ^ for plane waves of cos- and sin-like type are 
defined by the following equation: 



The spectral-density coefficients a k c (s) ^ n a multicomponent system, in contrast to 
a J k , instead of being normalized to unity are normalized by the following value: 

^ | - k ' c ( s ) 1 ~ V m-lw, ,J 2 ' [ J 

j=l iLi "M w k,c(s) | 

having the value ~ 0.7 in the case of vitreous silica. 

Eqs. (10) - (15) fully determine the time evolution of different k' plane-wave com- 
ponents in the propagating vibrational excitation via the spectral densities and the 
vibrational spectrum itself. Another useful characteristic often used to characterize 
the decay of an initial external excitation is the time correlation function [36], 

(u(t) -u(O)) _ EiU,(t)u;(0) _ EaK cos ^ ( 16 ) 



(u(O)-u(O)) " £iW kii w kii ^ (w 2 > 

where the spectral-density coefficient a k is defined by Eq. (14) with w k)C ( s ) replaced 
by w k , and is related to the spectral-density coefficient a{. according to the following 
equation: 

a{ = Y / a((e j m~ 1 ei'} . ( 17 ) 

r 

In the case of a one-component system, the spectral-density coefficients a k and a k 
are obviously identical. 

The decay of the external plane- wave excitation can also be characterized via the 
properties of the final state after decay, averaged over time as t — > oo. An initial 
plane-wave excitation characterized by the wavevector k and polarization n is scat- 
tered to different plane-wave components characterized by the wavevectors k' and 
polarizations n'. The distribution, p(k', n'|k, n), of the weights of different plane- 
wave components averaged over time in the final state, a k / fi /(t), is of particular 
interest. Such a distribution is defined by the following equation: 



1 ( Ej Kl 2 |q k ', s l 2 Ej \®U 2 hl 

2 (W> s )) 2 ' «w 2 » 2 



p(k', n'|k, n) = al(t) = ± ^ ' \v> + V,,. \v>' > ( 18 ) 



where, as in Eq. (1) and subsequently, the wavevector index k' also includes the 
polarization index n'. Eq. (18) is obtained by averaging Eq. (10) over time and 
using Eq. (12). Bearing in mind that (w k , s ) ~ (w k , c ) ~ (w k ,) and that the sum 

\<xh sP + l^k' c 1 2 i s independent of the phase of cos- and sin-like components defined 
by Eq. (9), expression (18) can be transformed to: 

V • W \ 2 \ry j I 2 

p(k-,n'|k,n) = ^^' . (19) 



The distribution (19) averaged over all polarizations n' in the final state is 

Pav(k'|k,n) ~2p(k / ,fi;|k,fi)+p(k / ,n;|k,n) , (20) 

where the unit vector nj. stands for transverse polarization in the final state while h[ 
refers to longitudinal polarization, and the factor 2 takes into account the existence 
of two independent and, in glasses, equivalent transverse polarizations. Glasses are 
isotropic, and an averaging of Eqs. (18)-(20) over the directions of both initial and 
final wavevectors (including averaging over transverse polarizations in the initial 
wave) can be made, resulting in: 

p w ,tQ)( k '\ k ) = (pav(k / |k,fi t(1) ))^ k , . (21) 

If we are interested in the contribution of the same plane-wave k-component as 
in the initial excitation, then the wavevector k' should be replaced by k in Eqs. 
(18)- (20) and averaging only over k-directions should be made in Eq. (21). 

Ill SPECTRAL DENSITIES 

As follows from Sec. II, the coefficients a J k , in the expansion of different k- 
plane-waves over the eigenmodes, i.e. projections of plane waves onto eigenvectors, 
and related spectral densities, \a^\ 2 , \al.\ 2 and a^a^, fully determine the dynamical 
response of the system to the initial plane-wave excitation. We calculate below the 
spectral densities for two models of disordered structures: (i) a model of v-SiC>2 
constructed by molecular dynamics and (ii) a disordered zig-zag chain. 

A Vitreous silica 

The models of v-SiO"2 have been constructed by iVPT-molecular-dynamics sim- 
ulations, using the potential of van Beest [37]. The van Beest potential has been 
modified for small interatomic distances according to Guissani and Guillot [38] . At 
large interatomic distances, we have used a cutoff for short-range interactions, mul- 
tiplying the modified van Beest potential by a Fermi-like step function. The step 
function is characterized by the step position at R cnt = 5.5A and the step width 
SRcut — 0.5A for all atomic species. The latter cutoff has been used to obtain a 
density of the glassy structure (at zero pressure), of 2.38g/cm 3 , reasonably close to 
the experimental value of 2.2g/cm 3 (see the discussion of the densification problem 
in Ref. [38]). Note that a similar cutoff (R cu t = 5.0 A and 5R cut = 0) has been used 
by Vollmayr et al. [39]. 

All glassy models have been created by quenching from the melt (T = 6000K) to 
the well-relaxed glassy state (T ~ 10~ 4 K) at an average quench rate of ~ lK/ps. 
No coordination defects have been found in the models. The fully dense dynamical 
matrices for the relaxed systems were diagonalized directly, resulting in eigenvec- 
tors {e J } and eigenvalues (oUj), thus allowing us to perform a complete harmonic 



vibrational analysis. Structural characteristics and vibrational properties of the 
models are very similar to those described in Ref . [40] . 

The models of v-SiC>2 were of two types: a cubic model containing N = 1650 
atoms and of box length L ~ 28. 4 A, and a bar configuration containing N = 1500 
atoms of size 85. 6A x 15. 6A x 15. 6A (a bar-shaped model of B2O3 has been also 
used in Ref. [41]). The bar-shaped models were constructed to allow access to much 
lower values of k (> 0.07A -1 ) for modes propagating along the bar than can be 
obtained for the cubic models (k > 0.22A -1 ). 

Silica is a multicomponent system so that the three spectral densities differ from 
each other due to different contributions in the coefficients and a£ of the mass 
factor (see Eqs. (3) and (14)). In the case of vitreous silica, the masses of the 
atomic species are quite comparable and the mass factor is of the order of unity, 
so that the different types of the spectral densities differ only slightly from each 
other and the following approximate relationships can be used: |«k| 2 — Ai|a{.| 2 and 
«k^k — ^2|«k| 2 ! "with the normalization constants for the corresponding spectral 
densities being A\ = J2j \oL\ 2 an d ^2 = Z)j«k^k> which, in the case of vitreous 
silica, give A 1 ~ 0.7 and A 2 ~ 0.8. 
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FIGURE 1. The spectral densities \a J k \ 2 for transverse ((a), (b) and (c)) and longitudinal ((d), 
(e) and (f)) initial polarizations for different magnitudes, k, of the initial wavevector as shown in 
the figure. 



The shape of the spectral density depends on the characteristics of the plane 
wave (wavevector and polarization) and on the atomic structure itself. In disor- 
dered structures and for small values of the wavevector magnitude (ka < 1), the 
spectral density both for longitudinal and transverse polarizations has the shape 
of a single pronounced peak (see Figs. 1(a) and 1(d)). The positions of these low- 
frequency peaks at v t ,\ are related to the wavevector magnitude according to the 
linear dispersion relation (see Fig. 2), 

i/ t> i ~ Ct,iA;/27r . (22) 

As seen from Fig. 2, the calculated dots in the low- frequency range lie on the straight 
lines plotted using experimental sound velocities (c t — 37.5A/ps and q ~ 59A/ps 
[42]). 

With an increase of the magnitude k of the wavevector, the peak-shaped spectral 
density shifts to higher frequencies and its width increases (see Figs. 1(b) and 1(e)). 
At large enough k > lA -1 , the spectral density no longer consists of a single peak 
but rather resembles the vibrational density of states (VDOS) (see Figs. 1(c) and 
1(f)), clearly showing the two frequency bands found in v-Si0 2 [40]. 
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FIGURE 2. The dispersion laws for transverse (a) and longitudinal (b) polarizations of the 
initial plane-wave excitation. The solid circles and squares were obtained from the fit of the 
spectral densities by Lorentzians and the DHO model, respectively. The stars in (b) correspond 
to IXS data by Benassi et al [42] . The solid lines represent the long- wavelength limit characterized 
by the experimental sound velocities. 



If the spectral densities are peak-shaped, two of their characteristics, the peak 
position and width, are normally used in order to describe the propagation of ex- 
ternal plane-wave excitations [19,21,42,43]. The peak position is associated with 
the average frequency of the propagating excitation, while the peak width is asso- 
ciated with the decay time of the excitation. Indeed, if we look at the k-plane-wave 
component in the propagating excitation, u(t), its evolution with time is described 
by relation (7) at k' = k. The weight (amplitude) of this component, a^(t), de- 
cays with time according to Eqs. (10) and (12). A rough estimate of the time 
dependence of a^it) can be obtained if we assume that a^it) oc a^ c (t), i.e. a^it) is 
approximately the back cosine Fourier transformation of the spectral density 
(see Eq. (12)). If the spectral density has the shape of a well-defined peak which 
can be fitted, say, by a Lorentzian, i.e. 

/L "-(c-^) 2 + (r./2) 2 ' {26) 

where the Lorentzian position and full-width at half-maximum (FWHM) T w are 
the fitting parameters, then the back cosine Fourier transform of the function (23) 
is 



a k (t) ~ A 2 coscj fc texp{-r w t/2} , (24) 

where we have actually used Eq. (23) to fit the spectral density |«kl 2 normalized to 
unity and then took into account the factor A 2 (see the beginning of the section). 
As clearly seen from Eq. (24), the decay of the k-plane-wave component can be 
characterized by the average radial frequency and the inverse decay time 

t-,7 1 ~ r w (fc)/2 = Trr^fc) , (25) 

with r„(THz) = r w /27r. 

In Refs. [42,21], the damped harmonic oscillator (DHO) model has been used 
to fit spectral densities, which gives similar values for the average frequency and 
width, if (r^) 2 <C uo\. This inequality holds true in the region k < lA^ 1 and in 
particular in the IR regime around k ~ O.lA -1 (see below), where the spectral 
densities have a well-defined peak shape and fitting of the spectral densities by 
Lorentzian and/or DHO curves makes sense. 

We have used fits both by the Lorentzian and DHO models to obtain the av- 
erage frequency and decay time (not shown, see Ref. [46] for more detail) of the 
propagating plane-wave excitation as a function of the initial wavevector. The de- 
pendence of = a7fc/27r vs. k shown in Fig. 2 can be associated with some sort 
of "dispersion law". Of course, the propagating plane- wave excitation cannot be 
characterized by only one wavevector (and single frequency) and instead consists 
of a packet of plane waves (see Eq. (6)) with different wavevectors (packet of eigen- 
modes characterized by different frequencies). We chose from the k'-packet only 
one component characterized by the same wavevector as that of the initial plane 



wave and followed its time evolution. In that case, the dependencies presented 
in Fig. 2 can be regarded as the dispersion laws for a single plane- wave component. 
The experimental data for longitudinal external plane-wave excitations from IXS 
experiments [42], obtained by fitting the experimental curves with the DHO model, 
are shown by the stars in Fig. 2(b) and they agree well with our results. 

Note that the dispersion laws for both branches are practically linear in the low- 
frequency (long-wavelength) regime for v < 3THz. Above this frequency, a sort of 
"fast-sound" behaviour is observed. The increase in the slope of Vk is related to the 
changes in the shape of the spectral densities. A shoulder on the high-frequency 
side of the spectral-density peak for the longitudinal branch starts to appear at 
k > 0.3A -1 {y > 0.3THz - see Figs. 1(d), (e)). A similar transformation happens 
with the peak for the transverse branch at k > 0.5A" 1 . 



B Disordered zig-zag chain 



The other structural model we consider here is a toy model, namely a linear zig- 
zag chain on the plane. Toy models are very useful in studying atomic dynamics. 
Usually, scalar models are investigated because of their simplicity [34,44,45]. How- 
ever, important effects related to mixing of modes of different polarizations [46] are 
missed in such models. That is why we consider below one of the simplest vectoral 
models, namely a zig-zag chain in the x — y plane (see Fig. 3). 

Consider a chain with 2 atoms in the unit cell: {xi,yx} and {x 2 , y 2 }, 
so that the positions of the other atoms can be found as rj ^ = + za and 

r i°2 = r 2 + ^ a w hh i — 0, ±1, . . ., where a = {a, 0} is a unit-cell vector, so that 
the chain is directed along the x axis (see Fig. 3). 
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FIGURE 3. Zig-zag chain with two atoms in the unit cell. 



The nearest neighbours of different types (atoms 1 and atoms 2) are connected 
by springs with spring constants k^J in the same unit cell i and with k^J in different 
unit cells (i and i + while the nearest neighbours of the same type are linked by 
springs with constants k^J for atoms of the first type (the lower row) and k^J for 
the atoms of the second type (upper row). 

The potential energy then has the following expression: 

nri,i.n,2,-) = ^E{«15(l^- r Ml-|rS-'ai) 2 

/ i 

+ k (1) (\r v I lr (0) r (0) \\ 2 

+ K 2,i [\ r i+l,l - r i,2\ - - r i,2 I ) 

+ E^O^-^-l-lr^-rg?!) 2 }, (26) 

3=1,2 

with i numbering all unit cells. 

Four branches occur in the dispersion relation: two acoustic (longitudinal and 
transverse) and two optic branches (see Fig. 4a). The analytical expression for 
the dispersion is available for the symmetric case (mi = m 2 = m; K^j = = 
K f} — K ¥J — K^;xi/a = yi/a = 0; x 2 /a = 0.5): 
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FIGURE 4. (a): Dispersion curves (u> = LO\Jmj 'kW vs k — k/(n/a)) for the zig-zag chain 
model characterized by the following parameters: the equilibrium coordinate of atoms in the 
unit cell, x\/a = 0; yi/a = 0, x 2 /a = 0.4; y 2 /a = 0.5, ratio of force constants k^/k^ = 0.1 
(k^ 1,2 ' = k^ 1 ' 2 ^), masses mi = m 2 = m and the total number of atoms TV = 2000. (b): The VDOS 
of the zig-zag chain model with the same set of parameters as in (a) (solid line) together with 
that for a disordered chain with fluctuations in force constants Sk^/k^ = 0.3, Sk^/k^ = 0.3 
(dashed line). 



^ = 1 + ^ ± [l + ^ ± A cos (ka/2)] 



(27) 



with A = 2 (1 — cos (ka)) (/c^/zcW). From this expression it is not difficult to 
obtain that in the low-frequency limit (uj — > 0) for the longitudinal acoustic branch 
not surprisingly u oc k, while for the transverse acoustic branch u oc k 2 . Such a 
/c 2 -dependence is typical for transverse vibrations of a linear chain and is related 
to the much smaller restoring force in the y-direction (compared to that in the 
x-direction) for long-wavelength vibrations because of the absence of the spring 
continuum in that direction. 

The vibrational density of states g(uS) contains two bands characterized by typ- 
ical van Hove singularities around the band boundaries (see Fig. 4b). The k 2 - 
dependence of the transverse acoustic branch results in an uj~ 1 / 2 singularity of the 
VDOS as the frequency approaches zero, u — > 0. 

We are interested mainly in disordered structures. Several different types of 
disorder can be introduced in the system: 

(i) Positional disorder which is due to random positional vectors rf\ 

(ii) Force-constant disorder: K^i £1X6 randomly distributed around their mean 
value, k, 

Ki = n + 5k- <5 iiran , (28) 

where 8k is a typical width of the distribution and 5i iran are random numbers 
distributed around zero with the density distribution, p(5), e.g. the normal distri- 
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bution, p(8) = exp {— 5 2 /2}/v / 27r , with unit variance. Possible negative values of 
the spring constants are replaced by their absolute values. 

(iii) Mass disorder: the atomic masses are randomly distributed, m ; = m + 5m ■ 

^j,rari' 

In the case of force-constant and mass disorder (which are considered in what 
follows), the equilibrium positions of the atoms are not changed. This is quite 
convenient and an ideal (crystalline) linear chain can be treated as the crystalline 
counterpart of the disordered chain. 

The VDOS of a disordered zig-zag chain is similar to that for the crystalline chain 
(see Fig. 4b) except all sharp features (excluding the region around zero) are washed 
out (band tails appear). We should also note that extra states (relative to the Debye 
spectrum) appear in the low- frequency regime (cu < c t) i(n/a), with c tj i being the 
transverse or longitudinal sound velocity. These extra states, characterized by the 
change in the VDOS, Ag(u) = g{uj) — g CTy st(u) , are called the boson peak [40]. 

The spectral densities for the disordered zig-zag chain, being of particular inter- 
est, are calculated according to Eq. (3) and presented in Fig. 5. 

As expected in the long- wavelength limit, the spectral density has the shape of 
a peak. The position of the peak can be approximately found from the dispersion 
for acoustic branches (see below). The width of the peak increases with increasing 
magnitude of the wavevector of the initial plane wave. This corresponds to more 
intense scattering of the plane waves by disorder. 



IV ANALYSIS OF THE FINAL STATE IN 
MOMENTUM SPACE 

As is well known, a scattering process can be investigated by analysis of the 
final state. First, we consider the final state of a single k-plane-wave component 
characterized by the same wavevector as the initial one. The phase of this wave has 
a random value and is not an informative characteristic. The important quantity 
is the amplitude of the wave, or more precisely its squared average value, defined 
by Eq. (18) with k' = k, 

- l [ E,I^I 2 K S | 2 Ej l4lVk,cl 2 l Ej \al\ 2 \al 
° k 2\ («,»' + (K iC »2 J" « W £»2 

This value can be easily estimated for a peak-shaped spectral density of width T. 
Indeed, the number of eigenmodes contributing to an initial plane wave is of order 
3N • (T/D), where D is the width of the whole vibrational spectrum (~ 40THz in 
the case of vitreous silica). Then we can easily evaluate from the normalization 
conditions Eqs. (5), (17) for the spectral densities the average value of the spectral 
density in the peak region, la 7 ^ 2 ~ |«k, s | 2 ~ {D/T) • (1/3 AT), and obtain the 
following estimate for a£, 



(29) 



where we have taken into account that (w 2 .) ~ 1 according to Eq. (5). The factor 
D/Y in relation (30) shows that the averaged squared amplitude is inversely pro- 
portional to the number of initially excited modes and not to all the modes. The 
factor D/Y in the Ioffe-Regel (IR) region is much larger than unity; D/Y ~ 10 2 
for T ~ ^ir/V ~ 0.3THz. Therefore we can say that the k-plane-wave component 
around and below the IR regime is not damped (the squared average amplitude is 
not of order 1/3N) but rather is attenuated (scattered), weakly (strongly) below 
(beyond) the IR limit as discussed below. The k-plane-wave component is damped 
at k > k* ~ lA^ 1 when the peak width becomes comparable to the full spectral 
width, r ~ D. 

The analysis of all k'-plane-wave components in the final state, in particular the 
distribution (18) of their weights, allows the scattering mechanism to be clarified. 
Let us consider an initial plane wave characterized by the wavevector k and polar- 
ization n. This wave is scattered with time into different plane waves characterized 
by wavevectors k' and polarizations n', which do not necessarily coincide with the 
initial polarization. We would like to know the weights of all plane-wave compo- 
nents in the final state as a function of the wavevector magnitude k'. The distribu- 
tions of the transverse and longitudinal plane waves, p(k', n(|k, n) and p(k', n(|k, n) 
(see Eq. (19)), and the distribution averaged over polarizations, p av (k'|k, n) (see 
Eq. (20)), for both transverse n t and longitudinal fii polarizations of the initial 
plane- wave excitation are of particular interest. These distributions depend only 
on the spectral densities |a^.| 2 , \q^,\ 2 and the vibrational spectrum itself, and can 
be easily calculated numerically for different k. 

A Vitreous silica 

The results of such calculations for vitreous silica are presented in Fig. 6. The 
upper (lower) row describes the scattering of initial longitudinal (transverse) plane 
waves, characterized by different wavevector magnitudes, into transverse and lon- 
gitudinal plane waves and also the distribution of the weights averaged over the 
polarization in the final state. 

First, we consider scattering of a longitudinal initial wave (the upper row in 
Fig. 6). The weight distributions, p(k', n(|k, fii) and p(k', n[|k, r^), characterize the 
scattering of the longitudinal wave to a longitudinal wave, the {/ — > 1} channel, and 
of the longitudinal wave to a transverse wave, the {/ — > t} channel, respectively. 
As follows from Fig. 6, these distributions are peak shaped but the positions of 
the peaks are different. The distribution for the {I — > 1} channel has a maximum 
around k' n ~ k\ = k (or maybe a bit below the initial wavevector), while the 
distribution for the {/ — > t} channel is mainly concentrated at a higher wavevector 
value, k[ t > k\. The distribution, p av (k'|k, fii), averaged over polarizations of the 
final state is a sum of double the distribution for the {/ — > t} channel plus the 



distribution for the {I — > 1} channel. If the peaks related to the individual channels 
and constituting the average distribution are narrow enough, then the distribution 
function p av (k'|k, fii) is doubly peaked (not clearly seen in Fig. 6). If the peaks 
are too wide, then p av (k'|k, ni) looks like a single wide peak (see Fig. 6) with a 
maximum position k[ m close to k[ t . 

Such a shape of the distributions of the weights of plane waves in the final state 
can be qualitatively understood in the following way. The distribution function 
p(k', n(|k, ni) of the transverse waves is an integral (sum in the case of a finite-size 
model) of the product of two spectral densities, la^nj 2 f° r longitudinal and fi , | 2 
for transverse polarization. Around the IR region and below it, these peak-shaped 
spectral densities have maxima at v\ ~ c\k/2ir and v[ ~ (\k! /2ir, respectively, which 
generally do not coincide with each other. Therefore, the distribution p(k', n(|k, fii) 
has a maximum around k[ t satisfying the equation, v\ ~ c\k/2n ~ c t k{ t /2iT ~ v' t , 
i.e. 

k' lt ~ cik/ct , (31) 

which is obviously greater than the wavevector of the initial longitudinal wave. 

The distribution of longitudinal waves for the {I — > 1} channel can be analysed 
in a similar manner. The main difference from the {I — > t} channel is that the 
spectral density of the longitudinal plane wave in the final state coincides with the 
spectral density of the initial longitudinal plane wave at approximately the same 
wavevector magnitude as for the initial wave, 
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FIGURE 6. The distribution functions p(k',n£|k, n) (circles), p(k', n(|k, n) (pluses) and 
p av (k'|k, n) (stars) for longitudinal ((a), (c) and (e)) and transverse ((b), (d) and (f)) initial 
polarizations of plane waves characterized by different initial wavevector magnitudes k for a 
structural model of v-SiC>2. 
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(32) 



Actually, the value k' n should be slightly shifted to lower values, because the height 
of the peak for the spectral density | 2 increases with decreasing k' and the 
maximum of the product of the spectral densities is reached in the low-frequency 
tail of the spectral density for the initial plane wave. 

The scattering of an initially transverse plane wave occurs similarly. In particular, 
the conclusion that the average frequency, i/, of the majority of the plane- wave 
components comprising the final state coincides with the average frequency, v, of 
the initial plane wave, 



holds true independently of the polarization of the initial plane-wave excitation. 
Therefore, we can roughly say that the disorder-induced scattering of the plane 
wave is approximately "elastic" (on average). This is not an absolutely precise 
conclusion because, first, the plane-wave components are distributed in frequency 
(composed of eigenmodes having different frequencies) in the initial and final states 
and, second, even the maximum of the distribution in the final state is slightly 
shifted to lower frequencies as compared to the initial one, as discussed above. 

In the case of the scattering of an initial transverse plane wave, two channels 
are available: {t — > 1} and {t — > t}. The distribution functions, p(k', n(|k, n t ) and 
p(k', n(|k, n t ), of the weights of plane waves in the final state for these channels 
have peaks located around the following values: 



As follows from Eq. (34) and Figs. 6(b),(d),(f), the peak for longitudinal waves lies 
below the initial k, while for transverse waves the peak approximately coincides 
with k, being slightly shifted to smaller values for reasons similar to those discussed 
above for the {/ — > /} channel. 

The distribution functions shown in Fig. 6 were obtained for a bar-shaped struc- 
tural model of v-SiC>2. Such a model is effectively one-dimensional and has restric- 
tions for the available initial k and final k' vectors, which are mainly directed along 
the bar in the low- A; limit. This also restricts the number of scattering channels. 
In order to check the influence of the dimensionality of the model on the scattering 
of plane waves, we have performed a similar analysis for a cubic (3-dimensional) 
model of v-SiC>2 and have not found any influence of the dimensionality of the 
model for the available wavevector magnitudes k > 0.22A~ 1 (for the cubic model). 

Poor statistics in the long-wavelength limit (see Fig. 6) is a finite-size effect re- 
lated to the restricted number of the wave- vectors allowed by the periodic boundary 
conditions. An analytical extrapolation approach [46] can be used to overcome such 
a shortcoming. 



v' ~ v , 



(33) 



k 



:' tl ~ c t k/c\ and k' tt ~ k . 



(34) 



B Zig-zag chain 



Another possible way to overcome the disadvantages of finite-size 3-D nu- 
merical models is to analyse low-dimensional models. Much lower wavevectors 
k > fc ffi n = 2ir/N 1 / d a are available, for example, in one-dimensional (d = 1) models 
as compared to the 3-D case, and the acoustic spectrum appears to be be much 
more dense. In order to check and support the analytical and numerical approaches 
presented above for the 3-D case, we have performed numerical experiments for a 
disordered 1-D zig-zag chain (see Sec. Ill B) and calculated the distribution function 
p(k', n'|k, n) for it. 

Our main purpose here is to calculate the distribution function p(k', n'|k, n) char- 
acterizing the scattering of a plane- wave excitation. First, we have calculated this 
distribution function for the crystalline counterpart (5ni = 0) and not surprisingly 
we found for k < 71/ a only {t — > t} and {/ — > 1} channels (see the lower row in 
Fig. 7 marking the peaks at k' tt = k and k' n = k for the {t — > t} and {I — > 1} 
channels, respectively). Disorder changes the situation dramatically and gives rise 
to the occurrence of {t — ► /} and {/ — > t} channels (see the upper row in Fig. 7), in 
complete agreement with the results of the fc-analysis given in Fig. 6 for the case 
of v-SiC>2. The positions of the additional peaks, at k' tl ({£ — > /} channel), and 
Kt ({^ ~~ > channel) can be obtained from the dispersion laws for the crystalline 
zig-zag chain by solving the equations: 
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FIGURE 7. The distribution function p(k', n'|k, n) for different scattering channels (as marked 
in the figure) for a disordered ((a) and (b)), Sn^ — Sk^ 2 ' — 1, 8m\ — <5m 2 = 0) and an ordered ((c) 
and (d)) zig-zag chain characterized by the following parameters: mi = 7712 = 1, K W = K ( 2) = 1, 
Xi/a = yi/a = 0, a^/a = 1/2/0- = 0.5. The initial wavevector magnitude ka/ir — 0.25. 



(35) 



and 

u\(k) = cu t (k' lt ) , (36) 

respectively (see arrows in Fig. 8). The width of the peaks increases with increasing 
disorder. We have also found a similar shape of the distribution function p (for four 
channels) for all wavevectors k < it /a with the corresponding u t and uj\ lying in 
the range of the dense spectrum. 



V SCATTERING MECHANISM 



From numerical calculations for both the 1-D zig-zag chain model and the 3- 
D model of v-Si0 2 , we have found that plane waves scatter not only to modes 
of approximately the same wavelength (and polarization) but also to modes of 
different wavelength (and polarization) but of similar frequency. The reason for 
such scattering is a natural question. 

In our simulations on zig-zag chains, we have found the {t — > 1} and {I — > t} scat- 
tering channels even for models not showing an appreciable increase of the VDOS 
in the low-frequency regime. Hence, the appearance of the {t — > /} and {I — > t} 
channels should be explained in terms of existing transverse and longitudinal acous- 
tic waves. Indeed, in the crystal, transverse and longitudinal acoustic phonons are 
orthogonal to each other, and hence transverse (longitudinal) plane waves cannot 
be scattered into longitudinal (transverse) plane waves (as we checked numerically; 
see the lower row in Fig. 7). Disorder leads to changes in the interaction energy, 




with the result that an acoustic phonon with a particular energy can couple to 
other phonons with closely comparable energies, including phonons with different 
polarizations and wavevectors. Therefore, the resulting eigenmodes contain com- 
ponents of different polarization and different wavevectors. A plane wave is not 
an eigenmode in the disordered stucture and it is scattered into eigenmodes of ap- 
proximately the same energy as that of the plane wave. These eigenmodes contain 
both transverse and longitudinal components and therefore an original plane wave, 
independent of its polarization, is scattered into both transverese and longitudinal 
plane waves. This gives a qualitative explanation of the existence of {t — > 1} and 
{/ — > t} scattering channels. 

In the 3-D case, the situation can be more complicated. Apart from the scattering 
mechanism due to the disorder-induced mixing of transverse and longitudinal plane 
waves discussed above, extra states (comprising the Boson peak) relative to the 
Debye spectrum (e.g. optic modes pushed down by disorder (see [46,47])) could 
participate in the hybridization between plane waves with different polarization. 
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